TimeSer

Периодограммы

Реальные данные

Считываем данные, строим исходный ряд и его периодограмму.

ts_data <- read.csv("MRTSSM4453USN.csv", header = TRUE, as.is = FALSE, sep = ',')[,1:2]
ts_data[,1] <- c(0:359)
ts1 <- ts_data[1:336,]
ts <- ts(ts_data[,2], start = c(1992, 1), frequency = 12)
ts1 <- ts(ts1[,2], start = c(1992, 1), frequency = 12)
ts2 <- ts(ts1[c(1:54)], start = c(1992, 1), frequency = 12)
plot.ts(ts1, type = 'l')

spectrum(ts1, log='no', method='pgram', detrend=FALSE)

spectrum(ts1, log='no', method='pgram')

На периодограмме все частоты имеют значительный вклад, в том числе присутствует частота 1/2 (6/12), то есть, возможно, мы встретим “пилу”.

Белый шум

Также продемонструруем периодограмму белого шума.

data <- rnorm(1000, 0, 1)
plot.ts(data, type="l")

spectrum(data, log='no')

spectrum(data, log='no', taper=0)

Частотные фильтры

Bandpass

Построим ряд, его периодограмму и применим фильтр для подавления частот вне диапазона 0.2–0.45.

plot.ts(ts1, type = 'l')

spectrum(ts1, log='no', method='pgram')

ts_data[1:336,]

##     DATE MRTSSM4453USN
## 1      0          1509
## 2      1          1541
## 3      2          1597
## 4      3          1675
## 5      4          1822
## 6      5          1775
## 7      6          1912
## 8      7          1862
## 9      8          1770
## 10     9          1882
## 11    10          1831
## 12    11          2511
## 13    12          1614
## 14    13          1529
## 15    14          1678
## 16    15          1713
## 17    16          1796
## 18    17          1792
## 19    18          1950
## 20    19          1777
## 21    20          1707
## 22    21          1757
## 23    22          1782
## 24    23          2443
## 25    24          1548
## 26    25          1505
## 27    26          1714
## 28    27          1757
## 29    28          1830
## 30    29          1857
## 31    30          1981
## 32    31          1858
## 33    32          1823
## 34    33          1806
## 35    34          1845
## 36    35          2577
## 37    36          1555
## 38    37          1501
## 39    38          1725
## 40    39          1699
## 41    40          1807
## 42    41          1863
## 43    42          1886
## 44    43          1861
## 45    44          1845
## 46    45          1788
## 47    46          1879
## 48    47          2598
## 49    48          1679
## 50    49          1652
## 51    50          1837
## 52    51          1798
## 53    52          1957
## 54    53          1958
## 55    54          2034
## 56    55          2062
## 57    56          1781
## 58    57          1860
## 59    58          1992
## 60    59          2547
## 61    60          1706
## 62    61          1621
## 63    62          1853
## 64    63          1817
## 65    64          2060
## 66    65          2002
## 67    66          2098
## 68    67          2079
## 69    68          1892
## 70    69          2050
## 71    70          2082
## 72    71          2821
## 73    72          1846
## 74    73          1768
## 75    74          1894
## 76    75          1963
## 77    76          2140
## 78    77          2059
## 79    78          2209
## 80    79          2118
## 81    80          2031
## 82    81          2163
## 83    82          2154
## 84    83          3037
## 85    84          1866
## 86    85          1808
## 87    86          1986
## 88    87          2099
## 89    88          2210
## 90    89          2145
## 91    90          2339
## 92    91          2140
## 93    92          2126
## 94    93          2219
## 95    94          2273
## 96    95          3265
## 97    96          1920
## 98    97          1976
## 99    98          2190
## 100   99          2132
## 101  100          2357
## 102  101          2413
## 103  102          2463
## 104  103          2422
## 105  104          2358
## 106  105          2352
## 107  106          2549
## 108  107          3375
## 109  108          2109
## 110  109          2052
## 111  110          2327
## 112  111          2231
## 113  112          2470
## 114  113          2526
## 115  114          2483
## 116  115          2518
## 117  116          2316
## 118  117          2409
## 119  118          2638
## 120  119          3542
## 121  120          2114
## 122  121          2109
## 123  122          2366
## 124  123          2300
## 125  124          2569
## 126  125          2486
## 127  126          2568
## 128  127          2595
## 129  128          2297
## 130  129          2401
## 131  130          2601
## 132  131          3488
## 133  132          2121
## 134  133          2046
## 135  134          2273
## 136  135          2333
## 137  136          2576
## 138  137          2433
## 139  138          2611
## 140  139          2660
## 141  140          2461
## 142  141          2641
## 143  142          2660
## 144  143          3654
## 145  144          2293
## 146  145          2219
## 147  146          2398
## 148  147          2553
## 149  148          2685
## 150  149          2643
## 151  150          2867
## 152  151          2622
## 153  152          2618
## 154  153          2727
## 155  154          2763
## 156  155          3801
## 157  156          2219
## 158  157          2316
## 159  158          2530
## 160  159          2640
## 161  160          2709
## 162  161          2783
## 163  162          2924
## 164  163          2791
## 165  164          2784
## 166  165          2801
## 167  166          2933
## 168  167          4137
## 169  168          2424
## 170  169          2519
## 171  170          2753
## 172  171          2791
## 173  172          3017
## 174  173          3055
## 175  174          3117
## 176  175          3024
## 177  176          2997
## 178  177          2913
## 179  178          3137
## 180  179          4269
## 181  180          2569
## 182  181          2603
## 183  182          3005
## 184  183          2867
## 185  184          3262
## 186  185          3364
## 187  186          3322
## 188  187          3292
## 189  188          3057
## 190  189          3087
## 191  190          3297
## 192  191          4403
## 193  192          2675
## 194  193          2806
## 195  194          2989
## 196  195          2997
## 197  196          3420
## 198  197          3279
## 199  198          3517
## 200  199          3472
## 201  200          3151
## 202  201          3351
## 203  202          3386
## 204  203          4461
## 205  204          2913
## 206  205          2781
## 207  206          3024
## 208  207          3130
## 209  208          3467
## 210  209          3307
## 211  210          3555
## 212  211          3399
## 213  212          3263
## 214  213          3425
## 215  214          3356
## 216  215          4625
## 217  216          2878
## 218  217          2916
## 219  218          3214
## 220  219          3310
## 221  220          3467
## 222  221          3438
## 223  222          3657
## 224  223          3454
## 225  224          3365
## 226  225          3497
## 227  226          3524
## 228  227          4681
## 229  228          2888
## 230  229          2984
## 231  230          3249
## 232  231          3363
## 233  232          3471
## 234  233          3551
## 235  234          3740
## 236  235          3576
## 237  236          3517
## 238  237          3515
## 239  238          3646
## 240  239          4892
## 241  240          2995
## 242  241          3202
## 243  242          3550
## 244  243          3409
## 245  244          3786
## 246  245          3816
## 247  246          3733
## 248  247          3752
## 249  248          3503
## 250  249          3626
## 251  250          3869
## 252  251          5124
## 253  252          3143
## 254  253          3212
## 255  254          3603
## 256  255          3464
## 257  256          3916
## 258  257          3776
## 259  258          3994
## 260  259          4056
## 261  260          3588
## 262  261          3741
## 263  262          4007
## 264  263          5147
## 265  264          3333
## 266  265          3261
## 267  266          3596
## 268  267          3643
## 269  268          4096
## 270  269          3966
## 271  270          4166
## 272  271          4139
## 273  272          3736
## 274  273          4003
## 275  274          4012
## 276  275          5444
## 277  276          3486
## 278  277          3397
## 279  278          3761
## 280  279          3768
## 281  280          4222
## 282  281          4104
## 283  282          4409
## 284  283          4140
## 285  284          3955
## 286  285          4145
## 287  286          4135
## 288  287          5634
## 289  288          3488
## 290  289          3642
## 291  290          3907
## 292  291          3966
## 293  292          4242
## 294  293          4307
## 295  294          4572
## 296  295          4307
## 297  296          4260
## 298  297          4261
## 299  298          4488
## 300  299          5812
## 301  300          3578
## 302  301          3606
## 303  302          4074
## 304  303          4077
## 305  304          4456
## 306  305          4482
## 307  306          4598
## 308  307          4452
## 309  308          4346
## 310  309          4343
## 311  310          4638
## 312  311          5972
## 313  312          3792
## 314  313          3792
## 315  314          4436
## 316  315          4143
## 317  316          4702
## 318  317          4740
## 319  318          4761
## 320  319          4697
## 321  320          4416
## 322  321          4555
## 323  322          4926
## 324  323          6128
## 325  324          3933
## 326  325          3916
## 327  326          4445
## 328  327          4358
## 329  328          4861
## 330  329          4769
## 331  330          4993
## 332  331          5017
## 333  332          4454
## 334  333          4676
## 335  334          5057
## 336  335          6326

bandpass_ts <- bandpass(ts_data[1:336,],fhigh=0.45,flow=0.05,win=2,p=.05,detrend = FALSE)

## 
## ----- BANDPASS FILTERING STRATIGRAPHIC SERIES-----
##  * Number of data points= 336 
##  * Sample interval= 1 
##  * Mean value removed= 3008.012

##  * Center of bandpass filter = 0.25 
##  * 269 pos/neg frequency pairs will be bandpassed

ts_b <- ts(bandpass_ts[,2], start = c(1992, 1), frequency = 12)
plot.ts(bandpass_ts[,2], type = 'l')
spectrum(ts_b, log='no', method='pgram')

Полученный ряд не имеет линейного тренда, так как тренд соответствует низким частотам, а также отсутствует гармоника с частотой 1/2 (\((-1)^n\)).

АЧХ

Построим АЧХ для нескольких фильтров:

afc <- function(filter, omega) {
k <- seq_along(filter) - 1
h <- function(o) sum(rev(filter) * exp(-k*1i * o)) 
abs(sapply(omega, h))
}

freq <- seq(0, pi, 0.001)
filt <- rep(1, 36)
omega <- freq/2/pi
plot(afc(filt, freq) ~ omega, type = 'l')

filt <- c(1,-1)
plot(afc(filt, freq) ~ omega, type = 'l')

Первый для скользящего среднего с окном 36, второй– переход к разностям.

Смоделируем сильно зашумленный ряд и попробуем выделить необходимую частоту.

Спектральная плотность

Применим скользящее среднее к периодограмме шума, для оценки спектральной плотности:

n <- 200
N <- 1000
wnoise <- rnorm(N, 0, 1)
plot.ts(wnoise, type="l")

B <- spectrum(wnoise, log='no', taper=0, plot=FALSE)
plot(B$freq, B$spec, type='l')
lines(B$freq, movavg(B$spec, n, type='t'), type = 'l', col = 'red')

#spectrum(wnoise, kernel = kernel("daniell", c(11,7,3)), log = "no")
#bandpass_ts2 <- bandpass(cbind(c(0:N-1), wnoise),flow=0.2,fhigh=0.5,win=2,p=.5,detrend = TRUE)
#ts_b2 <- ts(bandpass_ts2[,2], start = c(1992, 1), frequency = 12)
#plot.ts(bandpass_ts2[,2], type = 'l')
#spectrum(ts_b2, log='no', method='pgram')

Для белого шума получается константа.

w0 <- wnoise[1]
wnoise <- wnoise[2:N]
cor <- 0.8
rnoise = Reduce(function(prev_v, next_v) cor * prev_v + next_v * sqrt(1 - cor^2), wnoise, w0, accumulate = T)

plot.ts(rnoise, type="l")

B <- spectrum(rnoise, log='no', taper=0, plot=FALSE)
plot(B$freq, B$spec, type='l')
lines(B$freq, movavg(B$spec, n, type='t'), type = 'l', col = 'red')

#bandpass_ts3 <- lowpass(cbind(c(0:N-1), rnoise),fcut=.05,win=2,p=.2,detrend = TRUE)
#ts_b3 <- ts(bandpass_ts3[,2], start = c(1992, 1), frequency = 12)
#plot.ts(bandpass_ts3[,2], type = 'l')
#spectrum(ts_b3, log='no', method='pgram')

Для красного шума убывающая кривая.

Выделение тренда

Скользящее среднее

Выделение тренда с помощью скользящего среднего с разной длиной окна (берем длину окна, кратную периоду ряда–12):

ma <- function(x, n){filter(x, rep(1 / n, n), sides = 2)}
plot.ts(ts1, type = 'l')
lines(x = time(ts1), y = ma(ts1, 24), type = 'l', col = 'blue')
lines(x = time(ts1), y = ma(ts1, 48), type = 'l', col = 'red')

Недостатком является необходимость достраивать края.

n <- 48
plot.ts(ts1, type = 'l')
lines(x = time(ts1), y = movavg(ts1, n, type='t'), type = 'l', col = 'red')
lines(x = time(ts1), y = ma(ts1, n), type = 'l', col = 'blue')

HP-фильтр

Применим HP-фильтр с параметром лямбда = 12960.

plot.ts(ts1, type = 'l')
hp_filt <- hpfilter(ts1, freq=12960, type="lambda")$trend
lines(hp_filt, type = 'l', col = 'red')

Результат похож на скользящее среднее, но тренд выделен и на краях. Недостаток- сложность подбора параметра.

Параметрические модели

time_ts1 <- time(ts1)
model <- nls(ts1 ~ b1*time_ts1+b2, start = list(b1 = 1, b2 = 3))
plot.ts(ts1, type = 'l')
lines(x = time_ts1, y = predict(model, ts1), type = 'l', col = 'red')

Приближение многочленом первой степени плохо описывает поведение на краях.

model <- nls(ts1 ~ b1*time_ts1^2+b2*time_ts1+b3, start = list(b1 = 1, b2 = 3, b3 = 0))
plot.ts(ts1, type = 'l')
lines(x = time_ts1, y = predict(model, ts1), type = 'l', col = 'red')

Приближение многочленом второй степени неплохо описывает тренд.

expR <- lm(log(ts1) ~ time(ts1))
plot.ts(ts1, type = 'l')
lines(x = time_ts1, y = exp(predict(expR)), type = 'l', col = 'red')

Приближение экспонентой неплохо описывает тренд.

SSA

Для выделения тренда рассмотрим метод ssa с длиной окна 12.

ssa.result <- ssa(ts1, L=12)
plot(ssa.result, type = 'vectors', idx=1:12)

plot(ssa.result, type = 'paired')

plot(ssa.result, type = 'series')

plot(ssa.result, type = 'wcor')

parestimate(ssa.result)

## $F1
##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf   0.003176 |  1.00318   0.00 |  1.00318   0.00000
## 
## $F2
##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf  -0.568060 |  0.56662   0.00 |  0.56662   0.00000
## 
## $F3
##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf  -0.757734 |  0.46873   0.00 |  0.46873   0.00000
## 
## $F4
##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf  -2.784229 |  0.06178   0.00 |  0.06178   0.00000
## 
## $F5
##    period     rate   |    Mod     Arg  |     Re        Im
##     2.000  -3.250759 |  0.03874   3.14 | -0.03874   0.00000
## 
## $F6
##    period     rate   |    Mod     Arg  |     Re        Im
##     2.000  -0.150974 |  0.85987   3.14 | -0.85987   0.00000
## 
## $F7
##    period     rate   |    Mod     Arg  |     Re        Im
##     2.000  -0.157958 |  0.85389   3.14 | -0.85389   0.00000
## 
## $F8
##    period     rate   |    Mod     Arg  |     Re        Im
##     2.000  -0.555808 |  0.57361   3.14 | -0.57361   0.00000
## 
## $F9
##    period     rate   |    Mod     Arg  |     Re        Im
##     2.000  -0.877564 |  0.41579   3.14 | -0.41579   0.00000
## 
## $F10
##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf  -0.113651 |  0.89257   0.00 |  0.89257   0.00000
## 
## $F11
##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf  -0.199074 |  0.81949   0.00 |  0.81949   0.00000
## 
## $F12
##    period     rate   |    Mod     Arg  |     Re        Im
##     2.000  -0.007688 |  0.99234   3.14 | -0.99234   0.00000

ssa.rec.trend <- reconstruct(ssa.result, groups = list(Trend = c(1)))
plot(ssa.rec.trend)

Трендом будем считать только первую компоненту.

plot.ts(ts1, type = 'l')
lines(ssa.rec.trend$Trend, col='red')

Из рассмотренных методов ssa наилучшим образом справляется с выделением тренда.

LOWESS

Добавим еще lowess

plot.ts(ts1, type = 'l')
#lines(lowess(ts1 ~ time(ts1), f=.05, iter=3L), col='green')
lines(lowess(ts1 ~ time(ts1), f=.1, iter=3), col='blue')
#lines(lowess(ts1 ~ time(ts1), f=.15, iter=3L), col=4)
lines(lowess(ts1 ~ time(ts1), f=.2, iter=3), col='red')

Результаты похожи на ssa, сильных отличий при использовании в lowess 10% и 20% точек нет.

Разложение на тренд, периодику и шум

Seasonal decomposition

dec_ts1 <- stats::decompose(ts1, type = 'multiplicative')
plot(ts1, type = 'l')
lines(dec_ts1$trend, col='red')

plot(dec_ts1$seasonal)

plot(dec_ts1$random)

spectrum(na.omit(dec_ts1$random), log='no', method='pgram', detrend=TRUE)

Полученный тренд неплохо описывает ряд, на периодограмме видно, что в шуме остались гармоники.

stl

Возьмем t.window=13, при общем количестве периодов 28:

stl_ts1 <- stl(log(ts1), s.window = 5, s.degree = 1, t.window = 13, t.jump = 1, t.degree=1)
spectrum(stl_ts1$time.series[,3], log='no', method='pgram', detrend=TRUE)

#acf(stl_ts1$time.series[,3])
#acf(stl_ts1$time.series[,2])
#acf(stl_ts1$time.series[,1])
plot(stl_ts1)

plot(ts1, type = 'l')
lines(exp(stl_ts1$time.series[,2]), col='red')

Полученный тренд неплохо описывает ряд, на периодограмме видно, что в шуме остались гармоники.

SSA

Рассмотрим ssa с максимально возможной длинной окна- половина ряда.

ssa.result <- ssa(ts1, L=336/2)
plot(ssa.result, type = 'vectors', idx = 1:20)

plot(ssa.result, type = 'paired', idx = 2:20)

plot(ssa.result, type = 'series', groups = 1:20)

plot(ssa.result, type = 'series', groups = list(1:20))

plot(ssa.result, type = 'wcor', groups = 1:20)

parestimate(ssa.result, group=list(2:12))

##    period     rate   |    Mod     Arg  |     Re        Im
##     2.999   0.002750 |  1.00275   2.10 | -0.50225   0.86790
##    -2.999   0.002750 |  1.00275  -2.10 | -0.50225  -0.86790
##     5.999   0.002663 |  1.00267   1.05 |  0.50123   0.86839
##    -5.999   0.002663 |  1.00267  -1.05 |  0.50123  -0.86839
##     2.400   0.002556 |  1.00256   2.62 | -0.86807   0.50157
##    -2.400   0.002556 |  1.00256  -2.62 | -0.86807  -0.50157
##     3.999   0.002536 |  1.00254   1.57 | -0.00055   1.00254
##    -3.999   0.002536 |  1.00254  -1.57 | -0.00055  -1.00254
##    11.968   0.002470 |  1.00247   0.52 |  0.86747   0.50245
##   -11.968   0.002470 |  1.00247  -0.52 |  0.86747  -0.50245
##    -2.000   0.002032 |  1.00203  -3.14 | -1.00203  -0.00000

ssa.rec.trend <- reconstruct(ssa.result, groups = list(Trend = c(1,13:15,18), Seas=c(2:12,16,17,19,20)))

Трендом считаем компоненты- 1, 13, 14, 15, 18; периодикой- 2-12, 16, 17, 19, 20; остальное- шум.

#plot(ssa.rec.trend$Seas)
plot.ts(ts1, type = 'l')
lines(ssa.rec.trend$Trend, col='red')

spectrum(residuals(ssa.rec.trend), log='no', method='pgram', detrend=TRUE)

acf(residuals(ssa.rec.trend))

Тренд выделен хорошо, остаток похож на красный шум, то есть в шуме не осталось тренда и/или периодик.

ЛРФ

Модельный ряд без шума

Модельный ряд, ранг 5 (2 корня от косинуса и 3 от квадратичного тренда).

a <- -1/600
b <- 12
c <- 5
n <- c(1:1000)
ts_mod <- ts(c*exp(a*n)*cos(2/b*pi*n) + (n/200)^2 + 100)
plot(ts_mod, type = 'l')

rank <- 5
ssa <- ssa(ts_mod, L = rank+1, svd.method = 'svd')
r0 <- reconstruct(ssa, groups = list(signal = 1:rank))$signal
plot(ssa)

n_eig <- list(1:rank)
l <- lrr(ssa, group = n_eig)
r1 <- Rssa::roots(l)
r1

## [1] 1.0000252+0.0000429i 1.0000252-0.0000429i 0.9999498+0.0000000i
## [4] 0.8645832+0.4991674i 0.8645832-0.4991674i

plot(l)

parestimate(ssa, group = n_eig)

##    period     rate   |    Mod     Arg  |     Re        Im
##  145168.017   0.000025 |  1.00003   0.00 |  1.00003   0.00004
## -145168.017   0.000025 |  1.00003  -0.00 |  1.00003  -0.00004
##       Inf  -0.000051 |  0.99995   0.00 |  0.99995   0.00000
##    12.000  -0.001667 |  0.99833   0.52 |  0.86458   0.49917
##   -12.000  -0.001667 |  0.99833  -0.52 |  0.86458  -0.49917

par <- parestimate(ssa, groups = n_eig, method = "esprit")

Полученные 5 корней изобразили на единичной окружности (1 кратности 3, отвечает за квадратичный тренд).

Решая систему получим коэфициенты ЛРФ и нарисуем ряд на фоне модельного.

p <- (par$moduli[1]+par$moduli[2]+par$moduli[3])/3
a <- p
ro <- par$moduli[4]
w <- 2*pi/(par$periods[4])
param <- matrix(c(0, 0, 1, 1, 0, 
                  a, a, a, ro*cos(w), ro*sin(w),
                  2^2*a^2, 2*a^2, a^2, ro^2*cos(2*w), ro^2*sin(2*w),
                  3^2*a^3, 3*a^3, a^3, ro^3*cos(3*w), ro^3*sin(3*w),
                  4^2*a^4, 4*a^4, a^4, ro^4*cos(4*w), ro^4*sin(4*w))
                  , 5, 5)
c <- solve(t(param), ts_mod[1:5])
c

## [1]  2.499999e-05  4.370836e-05  1.000000e+02  4.322916e+00 -2.495837e+00

t <- ts(Re((c[3]+c[2]*n+c[1]*n^2)*a^n+c[4]*ro^n*cos(w*n)+c[5]*ro^n*sin(w*n)))
plot(t)
lines(ts_mod, col = 'green')

Полученный ряд хорошо описывает исходный.

Модельный ряд с шумом

Модельный ряд с шумом, ранг 5 (2 корня от косинуса и 3 от квадратичного тренда).

ts_mod <- ts_mod + rnorm(1000, 0, 0.15)
plot(ts_mod)

ssa_signal <- ssa(ts_mod, svd.method = 'svd')
r0 <- reconstruct(ssa_signal, groups = list(signal = 1:rank))$signal
plot(ssa_signal, type = 'vectors')

В соответсвии с результатми ssa, ранг = 4, видимо задуманный квадратичный тренд был описан комбинацией двух экспонент.

rank <- 4
n_eig <- list(1:rank)
l1<- lrr(ssa_signal, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]

l2 <- lrr(ssa_signal, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
par <- parestimate(ssa_signal, groups = n_eig, method = "esprit")

Отобрали 4 сигнальных корня.

par

##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf   0.000623 |  1.00062   0.00 |  1.00062   0.00000
##       Inf  -0.000821 |  0.99918   0.00 |  0.99918   0.00000
##    12.000  -0.001682 |  0.99832   0.52 |  0.86457   0.49916
##   -12.000  -0.001682 |  0.99832  -0.52 |  0.86457  -0.49916

o <- order(abs(par$periods), decreasing = TRUE)
periods <- (par$periods[o])
moduli <- par$moduli[o]
len <- length(ts_mod)
vars <- matrix(nrow = len, ncol = rank)
for (i in 1:rank) {
  if (abs(periods[i]) > 10000)
    vars[, i] <- moduli[i]^(1:len)
  else if (periods[i] == 2)
    vars[, i] <- (-moduli[i])^(1:len)
  else if (periods[i] > 0)
    vars[, i] <- 
      moduli[i]^(1:len) * sin(2 * pi * (1:len) / periods[i])
  else
    vars[, i] <- 
      moduli[i]^(1:len) * cos(2 * pi * (1:len) / periods[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
coefs0

##           X1           X2           X3           X4 
## 56.876944333 43.105170602  0.005296981  5.005089155

Нашли коэффициенты ЛРФ, как коэффициенты линейной регрессии.

Нарисуем полученный ряд на фоне модельного.

a1 <- par$moduli[1]
a2 <- par$moduli[2]
ro <- par$moduli[3]
w <- 2*pi/(par$periods[3])
t <- ts(coefs0[1]*a1^n+coefs0[2]*a2^n+coefs0[4]*ro^n*cos(w*n)+coefs0[3]*ro^n*sin(w*n))
plot(t)
lines(ts_mod, col = 'green')

Полученный ряд хорошо описывает исходный.

Реальные данные

Применим этот подход к реальному ряду.

plot.ts(ts1, type = 'l')

spectrum(ts1, log='no', method='pgram')

ssa.result <- ssa(ts1, L=336/2)
plot(ssa.result, type = 'vectors', idx = 1:20)

plot(ssa.result, type = 'paired', idx = 2:20)

plot(ssa.result, type = 'series', groups = 1:20)

plot(ssa.result, type = 'series', groups = list(1:20))

plot(ssa.result, type = 'wcor', groups = 1:20)

Ранг равен 15, тренд описывается комбинацией линейной функции и гармоники с большим периодом, сезонность описывается 11 компонентами (пары периодов + “пила”).

rank <- 15
r0 <- reconstruct(ssa.result, groups = list(signal = 1:rank))$signal
n_eig <- list(1:rank)
l1<- lrr(ssa.result, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]

l2 <- lrr(ssa.result, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
main.roots

##  [1]  1.0029434+0.0000000i -0.5021671+0.8679024i -0.5021671-0.8679024i
##  [4]  0.5012392+0.8684362i  0.5012392-0.8684362i  0.8676007+0.5026391i
##  [7]  0.8676007-0.5026391i -0.0006864+1.0026209i -0.0006864-1.0026209i
## [10] -0.8678587+0.5016418i -0.8678587-0.5016418i -1.0016702+0.0000000i

plot(main.roots)

par <- parestimate(ssa.result, groups = n_eig, method = "esprit")
par

##    period     rate   |    Mod     Arg  |     Re        Im
##  6333.541   0.004064 |  1.00407   0.00 |  1.00407   0.00100
## -6333.541   0.004064 |  1.00407  -0.00 |  1.00407  -0.00100
##     2.999   0.002845 |  1.00285   2.10 | -0.50232   0.86797
##    -2.999   0.002845 |  1.00285  -2.10 | -0.50232  -0.86797
##     5.999   0.002752 |  1.00276   1.05 |  0.50124   0.86849
##    -5.999   0.002752 |  1.00276  -1.05 |  0.50124  -0.86849
##    11.965   0.002649 |  1.00265   0.53 |  0.86756   0.50264
##   -11.965   0.002649 |  1.00265  -0.53 |  0.86756  -0.50264
##     2.400   0.002638 |  1.00264   2.62 | -0.86815   0.50160
##    -2.400   0.002638 |  1.00264  -2.62 | -0.86815  -0.50160
##     3.999   0.002626 |  1.00263   1.57 | -0.00058   1.00263
##    -3.999   0.002626 |  1.00263  -1.57 | -0.00058  -1.00263
##     2.000   0.002112 |  1.00211   3.14 | -1.00211   0.00000
##    92.142  -0.002090 |  0.99791   0.07 |  0.99559   0.06800
##   -92.142  -0.002090 |  0.99791  -0.07 |  0.99559  -0.06800

Отобрали сигнальные корни.

p <- (par$moduli[1]+par$moduli[2])/2
o <- order(abs(par$periods), decreasing = TRUE)
periods <- (par$periods[o])
moduli <- par$moduli[o]
len <- length(ts1)
vars <- matrix(nrow = len, ncol = rank)
n <- (1:len)
for (i in 1:rank) {
  if (abs(periods[i]) > 1000)
    vars[, i] <- n^(i-1)*p^n
  else if (periods[i] == 2)
    vars[, i] <- (-moduli[i])^n
  else if (periods[i] > 0)
    vars[, i] <- 
      moduli[i]^n * sin(2 * pi * n / periods[i])
  else
    vars[, i] <- 
      moduli[i]^n * cos(2 * pi * n / periods[i])
}
lm0 <- lm(r0[1:len] ~ 0 + ., data = data.frame(vars))
coefs0 <- coef(lm0)
coefs0

##          X1          X2          X3          X4          X5          X6 
## 1663.215990   -1.273498   59.028995   41.087106   69.307551  -84.709186 
##          X7          X8          X9         X10         X11         X12 
##  171.204057  -48.993611  155.382530  -39.076783    3.108672  126.328773 
##         X13         X14         X15 
##    9.718763  150.410704   57.620405

Нашли коэффициенты ЛРФ, как коэффициенты линейной регрессии.

Нарисуем полученный ряд на фоне восстановленного.

ro1 <- par$moduli[7]
w1 <- 2*pi/(par$periods[7])
ro2 <- par$moduli[5]
w2 <- 2*pi/(par$periods[5])
ro3 <- par$moduli[11]
w3 <- 2*pi/(par$periods[11])
ro4 <- par$moduli[3]
w4 <- 2*pi/(par$periods[3])
ro5 <- par$moduli[9]
w5 <- 2*pi/(par$periods[9])
ro6 <- par$moduli[13]
w6 <- 2*pi/(par$periods[13])
ro7 <- par$moduli[15]
w7 <- 2*pi/(par$periods[15])
t <- ts((coefs0[1]+coefs0[2]*n)*p^n+
        coefs0[3]*ro7^n*cos(w7*n)+coefs0[4]*ro7^n*sin(w7*n)+
        coefs0[5]*ro1^n*cos(w1*n)+coefs0[6]*ro1^n*sin(w1*n)+
        coefs0[7]*ro2^n*cos(w2*n)+coefs0[8]*ro2^n*sin(w2*n)+
        coefs0[9]*ro3^n*cos(w3*n)+coefs0[10]*ro3^n*sin(w3*n)+
        coefs0[11]*ro4^n*sin(w4*n)+coefs0[12]*ro4^n*cos(w4*n)+
        coefs0[13]*ro5^n*sin(w5*n)+coefs0[14]*ro5^n*cos(w5*n)+
        coefs0[15]*ro6^n*cos(w6*n)
        , start = c(1992, 1), frequency = 12)
plot(r0)
lines(t, col = 'green')

plot(ts1)
lines(t, col = 'red')

Полученный ряд хорошо описывает восстановленный ряд.

Автоматический выбор трендовых компонент

Для начала зададим число групп, первый вариант– 4 групп (1 отвечает за тренд и 3 за сезонность), второй вариант– 6 групп (1 отвечает за тренд и 5 за сезонность),

plot(ssa.result, type = "vectors", idx = 1:20)

plot(ssa.result, type = "series", groups = 1:20)

plot(ssa.result, type = "wcor", groups = 1:20)

g <- grouping.auto(ssa.result, freq.bins = 4, threshold = 0)
plot(reconstruct(ssa.result, groups = g))

g <- grouping.auto(ssa.result, freq.bins = 6, threshold = 0)
plot(reconstruct(ssa.result, groups = g))

В первом случае в тренд попала сезонность, значит количество групп надо увеличить, во втором случае тренд выделен хорошо.

Второй случай– зададим freq.bins = 1/(12*2) (меньше 1/12- самой маленькой возможной частоты из сезонности) и threshold = 0 для определения порога.

g <- grouping.auto(ssa.result, freq.bins = list(Trend = 1/(12*2)), threshold = 0)
plot(reconstruct(ssa.result, groups = g))

plot(g, type='b', order = TRUE)

Будем считать компоненты частью тренда при привышении порога в 0.85.

g <- grouping.auto(ssa.result, freq.bins = list(Trend = 1/(12*2)), threshold = 0.85)
plot(ts1)
lines(reconstruct(ssa.result, groups = g)$Trend, col='red')

g$Trend

## [1]  1 13 14 15 18 31 32

К тренду отнесены компоненты 1, 13, 14, 15, 18, 31, 32. Тренд выделен хорошо.

Улучшение разделимости

fossa

fossa помогает, если есть слабая разделимость, но нет сильной. Здесь смешаны компоненты 13-15.

f <- fossa(ssa.result, nested.group = 13:15)

plot(ssa.result, type = "vectors", idx = 1:15)

plot(f, type = "vectors", idx = 1:15)

plot(ssa.result, type = "series", groups = 1:15)

plot(f, type = "series", groups = 1:15)

plot(ssa.result, type = "wcor", groups = 1:15)

plot(f, type = "wcor", groups = 1:15)

Не помогло :(

iossa

iossa уже и с отсутствием слабой разделимости может справиться. Подаем на вход компоненты, отнесенные к тренду и компоненты, отнесенные к периодике.

g <- list(c(1, 13:15), 2:12)
issa <- iossa(ssa.result, nested.group = g)

plot(ssa.result, type = "vectors", idx = 1:15)

plot(issa, type = "vectors", idx = 1:15)

## Warning in .contribution(x, idx, ...): Elementary matrices are not F-orthogonal
## (max F-cor is -0.000776). Contributions can be irrelevant

plot(ssa.result, type = "series", groups = 1:15)

plot(issa, type = "series", groups = 1:15)

r0 <- reconstruct(issa, groups = issa$iossa.groups)
plot(r0)

Собственные вектора 13-15 сгладились.

eossa

eossa также справляется с отсутствием сильной разделимости, на вход подаем количество групп = 8.

e <- eossa(ssa.result, nested.group = g, k = 8)

plot(ssa.result, type = "vectors", idx = 1:15)

plot(e, type = "vectors", idx = 1:15)

## Warning in .contribution(x, idx, ...): Elementary matrices are not F-orthogonal
## (max F-cor is -0.0119). Contributions can be irrelevant

plot(e, type = "series", groups = e$iossa.groups)

r0 <- reconstruct(e, groups = e$iossa.groups)
plot(r0)

Собственные вектора 13-15 также сгладились, в группы были объедены компоненты, отвечающие за линейную составлющую тренда, за гармонику с большим периодом из тренда и все периодики из сезонности попарно.

Прогноз

Рассмотрим два варианта: “отрежем” один и два периода и будем прогнозировать по оставшимся точкам “отрезанные”.

ts1s <- ts(ts1[1:324], start = c(1992, 1), frequency = 12)
ts1s2 <- ts(ts1[1:312], start = c(1992, 1), frequency = 12)
ssa.result_1 <- ssa(ts1s, L=324/2)
ssa.result_2 <- ssa(ts1s2, L=312/2)
plot(ssa.result_1, type = 'vectors', idx = 1:20)

plot(ssa.result_1, type = 'paired', idx = 2:20)

plot(ssa.result_1, type = 'series', groups = 1:20)

plot(ssa.result_1, type = 'series', groups = list(1:20))

plot(ssa.result_1, type = 'wcor', groups = 1:20)

plot(ssa.result_2, type = 'vectors', idx = 1:20)

plot(ssa.result_2, type = 'paired', idx = 2:20)

plot(ssa.result_2, type = 'series', groups = 1:20)

plot(ssa.result_2, type = 'series', groups = list(1:20))

plot(ssa.result_2, type = 'wcor', groups = 1:20)

При “отрезании” одного периода при группировке будем рассматривать 15 компонент.

При “отрезании” двух периодов при группировке будем рассматривать 12 компонент, так как дальше компоненты сильно смешаны.

Рекуррентный прогноз:

r <- rforecast(ssa.result_1, groups = list(1:15), len = 12)
plot(ts(c(ts1s,r), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')

rec_1 <- rmse(ts1[325:336], r)
rec_1

## [1] 108.9153

r <- rforecast(ssa.result_2, groups = list(1:12), len = 24)
plot(ts(c(ts1s2,r), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')

rec_2 <- rmse(ts1[313:336], r)
rec_2

## [1] 134.0593

Ошибка рекуррентного прогноза для 1 периода– 108.92, для 2 периодов– 134.06.

Векторный прогноз

v <- vforecast(ssa.result_1, groups = list(1:15), len = 12)
plot(ts(c(ts1s,v), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')

vec_1 <- rmse(ts1[325:336], v)
vec_1

## [1] 85.11311

v <- vforecast(ssa.result_2, groups = list(1:12), len = 24)
plot(ts(c(ts1s2,v), start=start(ts1), frequency=frequency(ts1)))
lines(ts1, col = 'red')

vec_2 <- rmse(ts1[313:336], v)
vec_2

## [1] 134.5759

Ошибка векторного прогноза для 1 периода– 85.11, для 2 периодов– 134.58.

Bootstrap

1. Для одного периода:

bf_r <- forecast(ssa.result_1, groups = list(1:15), method = "recurrent", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')

brec_1 <- rmse(ts1[325:336], bf_r$mean)
brec_1

## [1] 108.9153

bf_v <- forecast(ssa.result_1, groups = list(1:15), method = "vector", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')

bvec_1 <- rmse(ts1[325:336], bf_v$mean)
bvec_1

## [1] 85.11311

Ошибка bootstrap прогноза для 1 периода:

Рекуррентный– 108.92.

Векторный– 85.11.

Обычные и рекуррентный и векторный прогнозы показали такие же результаты.

Также построены 95%-предсказательные интервалы для обоих прогнозов.

2. Для двух периодов:

bf_r <- forecast(ssa.result_2, groups = list(1:12), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')

brec_2 <- rmse(ts1[313:336], bf_r$mean)
brec_2

## [1] 134.0593

bf_v <- forecast(ssa.result_2, groups = list(1:12), method = "vector", bootstrap = TRUE, len = 24, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')

bvec_2 <- rmse(ts1[313:336], bf_v$mean)
bvec_2

## [1] 134.5759

Ошибка bootstrap прогноза для 2 периода:

Рекуррентный– 134.06.

Векторный– 134.58.

Нет особых различий между векторным и рекуррентным прогнозами.

Обычные и рекуррентный и векторный прогнозы показали такие же результаты.

Также построены 95%-предсказательные интервалы для обоих прогнозов.

iossa для прогноза

На примере для двух периода:

Плохо отделяется тренд, для исправления этого используем iossa.

g <- list(c(1, 13:15), 2:12)

issa <- iossa(ssa.result_2, nested.group = g)

plot(ssa.result_2, type = "vectors", idx = 1:15)

plot(issa, type = "vectors", idx = 1:15)

## Warning in .contribution(x, idx, ...): Elementary matrices are not F-orthogonal
## (max F-cor is -0.00102). Contributions can be irrelevant

plot(ssa.result_2, type = "series", groups = 1:15)

plot(issa, type = "series", groups = 1:15)

Первые 4 собственных вектора, соответствующих компонентам тренда стали более гладкие.

Сравним прогноз тренда до и после.

Для сравнения выделим тренд каким-то другим способом на полных данных (рассмотрим ранее полученные результаты при помощи ssa на полных данных).

До:

bf_r <- forecast(ssa.result_2, groups = list(c(1, 13:15)), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_r)
lines(ts1, col = 'black', t='l')
lines(ssa.rec.trend$Trend, col='red', t='l')

brec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_r$mean)
brec_2

## [1] 125.9559

bf_v <- forecast(ssa.result_2, groups = list(c(1, 13:15)), method = "vector", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_v)
lines(ts1, col = 'black', t='l')
lines(ssa.rec.trend$Trend, col='red')

bvec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_v$mean)
bvec_2

## [1] 80.70864

После:

bf_r <- forecast(issa, groups = list(1:4), method = "recurrent", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_r)
lines(ts1, col = 'black', t='l')
lines(ssa.rec.trend$Trend, col='red')

brec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_r$mean)
brec_2

## [1] 172.6782

bf_v <- forecast(issa, groups = list(1:4), method = "vector", bootstrap = TRUE, len = 24, R = 100, level = 0.95)
plot(bf_v)
lines(ssa.rec.trend$Trend, col='red')

bvec_2 <- rmse(ssa.rec.trend$Trend[313:336], bf_v$mean)
bvec_2

## [1] 133.5523

Стало хуже :)

На восстановленных рядах 5, 7, 11, 13 виден резкий скачок в конце, а полученный тренд как раз находится значительно ниже тренда, полученного ssa на полных данных. Возможно причина именно в этом.

Прогноз с точки зрения корней ЛРФ

1. Для 1 периода:

rank <- 15
r0 <- reconstruct(ssa.result_1, groups = list(signal = 1:rank))$signal
n_eig <- list(1:rank)
l1<- lrr(ssa.result_1, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]

l2 <- lrr(ssa.result_1, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
main.roots

##  [1] -0.5021028+0.8682899i -0.5021028-0.8682899i  1.0028899+0.0000000i
##  [4]  0.8676498+0.5026699i  0.8676498-0.5026699i  0.5011646+0.8684447i
##  [7]  0.5011646-0.8684447i -0.8681241+0.5014793i -0.8681241-0.5014793i
## [10] -0.0006135+1.0023935i -0.0006135-1.0023935i -1.0016530+0.0000000i

par <- parestimate(ssa.result_1, groups = n_eig, method = "esprit")
par

##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf   0.006012 |  1.00603  -0.00 |  1.00603  -0.00000
##       Inf   0.003966 |  1.00397   0.00 |  1.00397   0.00000
##     2.999   0.003149 |  1.00315   2.10 | -0.50229   0.86834
##    -2.999   0.003149 |  1.00315  -2.10 | -0.50229  -0.86834
##     2.400   0.002799 |  1.00280   2.62 | -0.86841   0.50147
##    -2.400   0.002799 |  1.00280  -2.62 | -0.86841  -0.50147
##    11.969   0.002789 |  1.00279   0.52 |  0.86776   0.50258
##   -11.969   0.002789 |  1.00279  -0.52 |  0.86776  -0.50258
##     5.999   0.002781 |  1.00278   1.05 |  0.50120   0.86855
##    -5.999   0.002781 |  1.00278  -1.05 |  0.50120  -0.86855
##     3.999   0.002534 |  1.00254   1.57 | -0.00048   1.00254
##    -3.999   0.002534 |  1.00254  -1.57 | -0.00048  -1.00254
##     2.000   0.002153 |  1.00216   3.14 | -1.00216   0.00000
##    92.145  -0.001857 |  0.99814   0.07 |  0.99582   0.06801
##   -92.145  -0.001857 |  0.99814  -0.07 |  0.99582  -0.06801

Тренд описывается четырьмя компонентами; корень, соответствующий линейному тренду имеет модуль больше 1, значит линейный тренд возрастает, а корень, соответствующий гармонике, относящейся к тренду имеет модуль меньше 1, значит колебания затухают. Все корни, соответствующие гармоникам имеют модули больше 1, а значит колебания усиливаются.

Проверим это, спрогнозировав на 20 периодов вперед:

bf_v <- forecast(ssa.result_1, groups = list(1:15), method = "vector", bootstrap = TRUE, len = 12*20, R = 100)
plot(bf_v)

Тренд возрастает, а колебания увеличиваются.

2. Для 2 периодов:

rank <- 12
r0 <- reconstruct(ssa.result_2, groups = list(signal = 1:rank))$signal
n_eig <- list(1:rank)
l1<- lrr(ssa.result_2, group = n_eig)
r1 <- Rssa::roots(l1)
main.roots <- r1[Mod(r1) > 0.999]

l2 <- lrr(ssa.result_2, reverse = TRUE, group = n_eig)
r2 <- Rssa::roots(l2)
main.roots <- c(main.roots, 1/r2[Mod(r2) > 0.999])
main.roots

##  [1]  1.0031701+0.0000000i -0.5017732+0.8682682i -0.5017732-0.8682682i
##  [4] -0.8684850+0.5012696i -0.8684850-0.5012696i  0.8675172+0.5026221i
##  [7]  0.8675172-0.5026221i  0.5012153+0.8682713i  0.5012153-0.8682713i
## [10] -0.0004278+1.0024724i -0.0004278-1.0024724i -1.0022579+0.0000000i

par <- parestimate(ssa.result_2, groups = n_eig, method = "esprit")
par

##    period     rate   |    Mod     Arg  |     Re        Im
##       Inf   0.003225 |  1.00323   0.00 |  1.00323   0.00000
##     2.999   0.003017 |  1.00302   2.09 | -0.50199   0.86837
##    -2.999   0.003017 |  1.00302  -2.09 | -0.50199  -0.86837
##     2.400   0.002904 |  1.00291   2.62 | -0.86864   0.50129
##    -2.400   0.002904 |  1.00291  -2.62 | -0.86864  -0.50129
##    11.971   0.002738 |  1.00274   0.52 |  0.86777   0.50246
##   -11.971   0.002738 |  1.00274  -0.52 |  0.86777  -0.50246
##     5.999   0.002707 |  1.00271   1.05 |  0.50125   0.86843
##    -5.999   0.002707 |  1.00271  -1.05 |  0.50125  -0.86843
##     2.000   0.002630 |  1.00263   3.14 | -1.00263   0.00000
##     3.999   0.002615 |  1.00262   1.57 | -0.00037   1.00262
##    -3.999   0.002615 |  1.00262  -1.57 | -0.00037  -1.00262

Тренд описывается одной компонентой, соответствующий корень имеет модуль больше 1, значит тренд возрастающий. Все корни, соответствующие гармоникам также имеют модули больше 1, а значит колебания усиливаются.

Проверим это, спрогнозировав на 20 периодов вперед:

bf_v <- forecast(ssa.result_2, groups = list(1:12), method = "vector", bootstrap = TRUE, len = 12*20, R = 100)
plot(bf_v)

Тренд возрастает, а колебания увеличиваются.

Модификации SSA

Toeplitz-ssa

Генерируем стационарный ряд + шум.

b <- 12
n <- 1:120
ts_mod <- ts(cos(2/b*pi*n) + 2 + rnorm(120,0,0.07))
plot(ts_mod)

ssa_toep1 <- ssa(ts_mod, L = 60, kind = '1d-ssa')
ssa_toep2 <- ssa(ts_mod, L = 60, kind = 'toeplitz-ssa')
plot(ssa_toep1, type = 'vectors', idx = 1:5)

plot(ssa_toep1, type = 'paired', idx = 2:5)

plot(ssa_toep1, type = 'series', groups = 1:5)

plot(ssa_toep1, type = 'series', groups = list(5))

plot(ssa_toep1, type = 'wcor', groups = 1:5)

plot(ssa_toep2, type = 'vectors', idx = 1:5)

plot(ssa_toep2, type = 'paired', idx = 2:5)

plot(ssa_toep2, type = 'series', groups = 1:5)

plot(ssa_toep2, type = 'series', groups = list(1:5))

plot(ssa_toep2, type = 'wcor', groups = 1:5)

Сигнальными компонентами считаем первые 3– первая компонента- тренд, вторая и третья- косинус.

ssa.rec_toep1 <- reconstruct(ssa_toep1, groups = list(Trend = 1, Seas=c(2:3)))
ssa.rec_toep2 <- reconstruct(ssa_toep2, groups = list(Trend = 1, Seas=c(2:3)))

Посмотрим насколько восстановленный ряд отличается от модельного ряда.

ssa_signal1 <- ssa.rec_toep1$Trend + ssa.rec_toep1$Seas
ssa_signal2 <- ssa.rec_toep2$Trend + ssa.rec_toep2$Seas
rmse(ts_mod, ssa_signal1)

## [1] 0.06722906

rmse(ts_mod, ssa_signal2)

## [1] 0.09872395

Basic-SSA показывает результаты лучше, чем Toeplitz-SSA, уже при таком небольшом уровне шума.

Двойное центрирование

ssa.res_2center <- ssa(ts1, L=336/2, row.projector = "center", column.projector = "center")
plot(ssa.res_2center, type = 'vectors', idx = 1:25)

plot(ssa.res_2center, type = 'paired', idx = 3:25)

plot(ssa.res_2center, type = 'series', groups = 1:25)

plot(ssa.res_2center, type = 'series', groups = list(1:25))

plot(ssa.res_2center, type = 'wcor', groups = 1:25)

Трендом считаем компоненты– 1, 2, 14, 15, 16, 19, 22; сезонностью– 3-13, 17, 18, 20, 21.

ssa.rec <- reconstruct(ssa.res_2center, groups = list(Trend = c(1,2,14:16,19,22), Seas=c(3:13,17:18,20:21)))
plot.ts(ts1, type = 'l')
lines(ssa.rec$Trend, col='red')
lines(ssa.rec.trend$Trend, col='green')

spectrum(residuals(ssa.rec), log='no', method='pgram', detrend=TRUE)

acf(residuals(ssa.rec))

Сравним тренд, выделеный классическим ssa и ssa с двойным центрированием: на концах заметны отличия, адекватнее тренд описывает классический ssa, возможно, это связано с тем, что данный ряд имеет нелинейный тренд, так как в нем присутствует гармоника с большим периодом.

Остаток похож на красный шум, но для этого пришлось брать большее число компонент, чем при классическом ssa.

Заполнение пропусков, разладка

gapfill

Добавим пропуски.

ts_gap <- ts1
gap <- 193:246
ts_gap[gap] <- NA
plot(ts_gap)

Рассмотрим длину окна такую что, вторая часть имеет не ноль векторов вложений.

ssa_gap_res <- ssa(ts_gap, L=72)
plot(ssa_gap_res, type = 'vectors', idx = 1:20)

plot(ssa_gap_res, type = 'paired', idx = 2:20)

plot(ssa_gap_res, type = 'series', groups = 1:20)

plot(ssa_gap_res, type = 'series', groups = list(1:20))

plot(ssa_gap_res, type = 'wcor', groups = 1:20)

Результат применения gapfill:

gr <- list(c(1:12, 13, 14, 19, 20))
g <- gapfill(ssa_gap_res, groups = gr, method = "sequential", 
             base = "reconstructed")
plot(ts_gap, col = "black")
lines(g, col = "red")

plot(ts1, col = "black")
lines(ts(g[gap], start = c(2008, 1), frequency = 12), col = "red")

rmse(g[gap], ts1[gap])

## [1] 99.42902

igapfill

ig <- igapfill(ssa_gap_res, groups = gr, 
               base = "reconstructed", maxiter = 2000)
plot(ts_gap, col = "black")
lines(ig, col = "red")

plot(ts1, col = "black")
lines(ts(ig[gap], start = c(2008, 1), frequency = 12), col = "red")

rmse(ig[gap], ts1[gap])

## [1] 93.25305

igapfill c maxiter=2000 показывает результат лучше (с точки зрения rmse), чем gapfill.

Cadzow

Аппроксимируем исходный ряд рядом конечного ранга:

cadz <- cadzow(ssa.result, rank = 15, tol = 1e-10)
plot(ts1)
lines(cadz, col = 'red')

rmse(ts1, cadz)

## [1] 79.63702

plot(ssa(cadz))

Получили ряд конечного ранга = 15.

Weights-Cadzow

Рассмотрим вариант с весами, для достижения большей точности.

L <- 336/2
K <- length(ts1) - L + 1
alpha <- 0.01
weights <- vector(len = K)
weights[1:K] <- alpha
weights[seq(K, 1, -L)] <- 1
ssa_cadz_weights <- ssa(ts1, L = L, column.oblique = "identity", 
           row.oblique = weights)
cadz_weights <- cadzow(ssa_cadz_weights, rank = 15, maxiter = 10)
plot(ts1)
lines(cadz_weights, col = 'red')

rmse(ts1, cadz_weights)

## [1] 68.53668

plot(ssa(cadz_weights))

rmse меньше по сравнению с обычным Cadzow.

Оптимальный подбор параметров

K.sliding <- 2
forecast.base.len <- 2*frequency(ts1s)
base.len <- length(ts1s)
sliding.len <- base.len - K.sliding - forecast.base.len + 1
ncomp <- 1:22
L.min <- 4*frequency(ts1s)
Ls <- seq(L.min, length(ts1s)/2, by = frequency(ts1s))
m0 <- forecast.sliding.rmse(ts1s,
                           K.sliding = K.sliding,
                           L = Ls, ncomp = ncomp,
                           method = "recurrent",
                           forecast.len = forecast.base.len,  
                           .progress = "none")
par <- optim.par(m0)
par$L.opt

## [1] 96

par$ncomp.opt

## [1] 19

par$m

##      
## L            1        2        3        4        5        6        7        8
##   48  299953.1 233286.9 209547.6 171339.2 147837.7 109897.4 104276.8 84959.12
##   60  300770.9 234748.7 209677.6 170746.5 147855.6 109503.2 104222.1 82925.58
##   72  302237.5 236795.8 210847.1 171382.3 149042.5 110470.1 105539.2 82910.69
##   84  304244.0 239164.3 212749.6 173051.5 150872.3 112173.2 107477.1 84272.42
##   96  306387.1 241653.4 214909.2 175252.9 153017.0 114315.6 109536.0 86265.41
##   108 307353.6 242671.5 215892.9 176053.8 153983.8 115229.7 110586.8 87142.54
##   120 307666.1 243108.7 216192.7 176229.1 154257.1 115393.7 111138.9 87432.56
##   132 308207.8 243738.7 216734.6 176719.8 154825.6 115917.6 111898.0 87566.09
##   144 309398.5 244830.0 218041.6 178012.8 156068.5 117283.2 113231.9 88410.93
##   156 309888.9 245423.7 218581.7 178471.3 156625.0 117575.2 113664.7 88719.02
##      
## L            9       10       11        12        13        14        15
##   48  49754.03 37248.92 18130.66 10231.960  9687.210  6494.669  6028.608
##   60  48578.32 33844.42 16511.54  9067.978  8032.647  5893.959  7281.058
##   72  49511.84 32891.48 17099.49  9871.432  7895.322  8241.512  8317.616
##   84  51214.10 33827.41 18672.66 11617.241  7961.980  9283.513  7855.896
##   96  53270.65 35385.83 20812.33 13823.281  8583.109 13106.236 12130.460
##   108 54502.91 36572.21 22129.72 15225.695  9593.236 17043.091 14907.630
##   120 55160.41 37004.88 22963.98 16007.474 12661.867 19583.756 26913.524
##   132 55860.01 37605.57 23905.51 16866.262 15618.963 25711.432 38682.329
##   144 56987.76 39293.25 25145.77 18147.784 19164.729 27267.362 31318.891
##   156 57173.76 39604.56 25442.52 18487.078 22806.867 26377.526 32160.950
##      
## L            16        17        18        19        20        21         22
##   48   9539.674  8417.259 10341.970 10267.722 10008.970 11075.830  12608.131
##   60   7585.594  6416.566  6662.913  7452.579  8225.342  7565.183  12687.272
##   72   6158.725  5421.843  5165.162  6876.015  6794.874  6636.124   6226.282
##   84   7319.978  6057.928  4984.602  4840.112  6632.341  5832.768   5442.083
##   96  11265.733  5398.528  4362.144  4176.981 11777.094 10947.351  10365.682
##   108 14219.807  9028.584  8169.678  6615.222  4930.718  5243.344   5939.900
##   120 21554.983 19104.077 17760.675 16135.923  4229.633  4930.668   8401.067
##   132 36407.412 29051.667 27696.440 27311.061 15029.839 12463.532  40713.008
##   144 28135.336 24524.378 22608.222 21664.438  4676.953 94202.460 173117.126
##   156 29603.955 24596.428 21956.432 21456.467  7286.060 67050.744 146725.683

matplot(Ls, sqrt(par$m), ylab = "", xlab = "Window lengths",
        type = "l", col = topo.colors(10))

Построим прогноз (векторный и рекуррентный bootstrap) с оптимальными параметрами:

ssa.result_optim <- ssa(ts1s, L=96)
plot(ssa.result_optim, type = 'vectors', idx = 1:19)

plot(ssa.result_optim, type = 'paired', idx = 2:19)

plot(ssa.result_optim, type = 'series', groups = 1:19)

plot(ssa.result_optim, type = 'series', groups = list(1:19))

plot(ssa.result_optim, type = 'wcor', groups = 1:19)

bf_r <- forecast(ssa.result_optim, groups = list(1:19), method = "recurrent", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_r)
lines(ts1, col = 'black', t='l')

brec_3 <- rmse(ts1[325:336], bf_r$mean)
brec_3

## [1] 74.5293

bf_v <- forecast(ssa.result_optim, groups = list(1:19), method = "vector", bootstrap = TRUE, len = 12, R = 100, level = 0.95, interval = 'prediction')
plot(bf_v)
lines(ts1, col = 'black', t='l')

bvec_3 <- rmse(ts1[325:336], bf_v$mean)
bvec_3

## [1] 72.34684

rmse векторного прогноза немного меньше.

Разладка

Разладка есть в данных (последние два периода).

plot(ts)

N <- length(ts)
ssa_full_res <- ssa(ts, L = N/2)
plot(ssa_full_res, type = 'vectors', idx = 1:20)

plot(ssa_full_res, type = 'paired', idx = 2:20)

plot(ssa_full_res, type = 'series', groups = 1:20)

plot(ssa_full_res, type = 'series', groups = list(1:20))

plot(ssa_full_res, type = 'wcor', groups = 1:20)

Рассмотрим первые 12 компонент:

rank <- 12
M <- 120
L <- M/4
hm <- hmatr(ts, B = M, T = 36, L = L, neig = rank)
plot(hm)

На матрице неоднородности видна разладка в последних периодах.